A new paradigm for generating high-quality cardiac pacemaker cells from mouse pluripotent stem cells

Cardiac biological pacing (BP) is one of the future directions for bradyarrhythmias intervention. Currently, cardiac pacemaker cells (PCs) used for cardiac BP are mainly derived from pluripotent stem cells (PSCs). However, the production of high-quality cardiac PCs from PSCs remains a challenge. Here, we developed a cardiac PC differentiation strategy by adopting dual PC markers and simulating the developmental route of PCs. First, two PC markers, Shox2 and Hcn4, were selected to establish Shox2:EGFP; Hcn4:mCherry mouse PSC reporter line. Then, by stepwise guiding naïve PSCs to cardiac PCs following naïve to formative pluripotency transition and manipulating signaling pathways during cardiac PCs differentiation, we designed the FSK method that increased the yield of SHOX2+; HCN4+ cells with typical PC characteristics, which was 12 and 42 folds higher than that of the embryoid body (EB) and the monolayer M10 methods respectively. In addition, the in vitro cardiac PCs differentiation trajectory was mapped by single-cell RNA sequencing (scRNA-seq), which resembled in vivo PCs development, and ZFP503 was verified as a key regulator of cardiac PCs differentiation. These PSC-derived cardiac PCs have the potential to drive advances in cardiac BP technology, help with the understanding of PCs (patho)physiology, and benefit drug discovery for PC-related diseases as well.


INTRODUCTION
Under physiological conditions, the cardiac conduction system (CCS) sends electrical impulses spontaneously and periodically to the atria and ventricles to trigger the heart to produce regular contractile activity. 1,2The CCS consists of the sinoatrial node (SAN), the atrioventricular node (AVN), the right and left bundle branches, and the Purkinje fiber network.The electrical impulses of the CCS originate from parenchymal cells known as PCs in SAN. 2,3A variety of factors such as myocardial ischemia, myocarditis, cardiomyopathy, hypertension, heart failure, gene mutation and cardiac degeneration can cause cardiac PC-related disorders, which manifest bradyarrhythmia, tachyarrhythmia, cardiogenic syncope, cardiogenic shock and even sudden cardiac death. 4At present, the treatment of severe cardiac PC-related disorders mainly relies on implantable electronic cardiac pacemakers (ECP), but ECP does not correct defects in the cardiac PCs themselves. 5Moreover, although great advances have been made since its first introduction to clinical practice in the 1950s, inevitably, ECP-related complications are still not uncommon, such as pneumothorax, cardiac perforation, and infection. 5Cardiac BP is one of the alternative strategies to arrhythmia treatment. 6urrently, cardiac BP can be achieved by gene-based, cell-based, hybrid gene-cell, and cellular reprogramming-based approaches. 7rdiac BP has been the focus and hot spot of cardiac pacemaking in recent years.
However, the development of cardiac BP technology was restricted by a number of factors, the most central of which was the availability of high-quality cardiac PCs.The number of cardiac PCs in SAN is extremely low, with only about 10,000 PCs in human SAN.In view of this, stem cell differentiation and somatic reprogramming are two alternative approaches to obtaining cardiac PCs, although both methods produce heterogeneous cell populations.To purify cardiac PCs from other cell types, either an antibody against a specific PC marker or a fluorescence reporter cell line indicating specific PC gene expression is needed.
By lineage tracing during embryo development, a set of genes, including Hcn4, Tbx3, Tbx5, Tbx18, Shox2, and Isl1, have been validated to relatively specifically express in the forming SAN either transiently or constitutively. 1,3Among these cardiac PCspecific expression genes, the short stature homeobox 2 (SHOX2) and the hyperpolarization and cyclic nucleotide 4 (HCN4) come to our attention.SHOX2 plays a pivotal role in SAN development and function. 8,9Mechanistically, SHOX2 prevents the formation of working myocardium by inhibiting Nkx2-5 expression while activating Hcn4, Isl1, and Tbx3 expression. 1,3Meanwhile, HCN4 is a widely accepted PC marker, which is responsible for generating funny current (I f ), the main driver of automatic electrical activity in cardiac PCs. 10 So far, four members of HCN have been identified, Hcn1-4, which express differentially in the heart, with Hcn4 highly expressed in the SAN followed by Hcn1. 11In Hcn4 knock-out (KO) mice, a 70% reduction in I f current was recorded with an overall 60% reduction of spontaneous heart rate. 12The properties of Shox2 and Hcn4 may be utilized to target cardiac PCs during differentiation.
PSCs, either derived from the epiblast of early embryos (embryonic stem cells, ESCs) or reprogrammed from somatic cells (induced pluripotent stem cells, iPSCs), can provide a cell source for cardiac BP by in vitro differentiation. 13Both cardiac PCs and working cardiomyocytes have been derived from PSCs.][16][17] Epiblast in blastocyst is the founder of the body, which generates three germ layers and germ cells.In parallel, naïve ESCs have been derived from epiblast in blastocyst (embryonic days 3.5-4.5,E3.5-E4.5)while primed PSCs, termed epiblast stem cells (EpiSCs) can be derived from post-implantation embryo epiblast (E5.5-E6.5). 18Unlike naïve ESCs, EpiSCs display restricted pluripotency and more heterogenous differentiation potential. 19,20t has been reported that in vitro, the differentiation efficiency from naïve ESCs directly is not as high as expected. 21Therefore, it has been postulated that there are intermediate stages between naïve and primed pluripotent states, and as naïve epiblast exits pluripotency, it undergoes a capacitation stage to acquire full differentiation potential, which has been corroborated by both in vivo analysis and in vitro derivation of formative stem cells (FSCs) from E5.5-E6.0post-implantation epiblast. 22Apart from the potential to generate primordial germ cells (PGCs), 23 FSCs can differentiate into cell types representing three germ layers with higher efficiency. 22,24ere, by applying PC-specific markers and simulating the developmental route of PCs, we developed an integrated strategy to generate cardiac PCs from PSCs.Shox2:EGFP; Hcn4:mCherry mESC reporter line was constructed to track the in vitro differentiation of cardiac PCs.By following the naïve to formative transition and manipulating the signaling pathways involved in cardiac PCs development, we established a simple and efficient cardiac PCs differentiation method to produce more than 20% SHOX2 + ; HCN4 + cardiac PCs.Both in vivo and in vitro characterization revealed the similarity between PSC-derived cardiac PCs and bona fide PCs.With scRNA-seq, we mapped the differentiation trajectory from naïve ESCs to cardiac PCs and validated Zfp503, a target of RA, 25 as a key factor regulating cardiac PCs differentiation.

RESULTS
Construction of Shox2:EGFP; Hcn4:mCherry mESC reporter line Shox2 and Hcn4 have been identified as two key genes contributing to cardiac PCs' gene regulatory network and electrophysiological properties, respectively. 26,27Overexpressing HCN4 in human iPSCs and mESCs-derived cardiomyocytes enables the cells to possess partial cardiac PC characteristics and exogenous SHOX2 expression promoted the differentiation of ESCs into cardiac PCs. 16,28,29To more precisely monitor the rising of cardiac PCs from PSCs and obtain functional PCs for cardiac BP, we set to construct an mESC reporter line with EGFP and mCherry fluorescence genes knocked into Shox2 and Hcn4 3' untranslated region (3'UTR) respectively by CRISPR/Cas9 mediated homologydirected repair (HDR) (Fig. 1a).In order to minimize side impact, we selected the protospacer-adjacent motif (PAM) on the 3'UTR of target genes close to the last exon to design single guide RNA (sgRNA) (Supplementary Fig. S1a).First, mCherry was introduced into 3'UTR of the Hcn4 locus in wild-type (WT) E14TG2a mESCs.
From the 50 clones we picked, by genomic PCR, 30 clones were identified as homozygous (HM) Hcn4:mCherry mESCs.We then selected and expanded one HM Hcn4:mCherry clone for Shox2:-EGFP targeting.We picked 64 clones for genotyping and identified five HM Shox2:EGFP; Hcn4:mCherry mESC clones (Supplementary Fig. S1b, c).Further Sanger sequencing of PCR products verified the in-frame sequence of the target gene and fluorescence reporter (Supplementary Fig. S1d).Furthermore, we analyzed some putative off-targets by PCR and sequencing, which revealed no mutations in these sites (Supplementary Fig. S1e, f).
While culturing the clones in naïve 2i/LIF medium, the colonies manifested dome-shaped morphologies, and no visible fluorescence was observed under a fluorescence microscope (Supplementary Fig. S1g).At the transcription level, most of the naïve pluripotency markers, including Oct4, Sox2, Nanog, and others, maintain their expression similar to those in WT mESCs (Fig. 1b).Immunostaining also confirmed homogenous expression of NANOG and OCT4 in Shox2:EGFP; Hcn4:mCherry mESCs (Fig. 1c).These results indicated that we have successfully established the Shox2:EGFP; Hcn4:mCherry mESC reporter line.
Validation of the credibility of the mESC reporter line for cardiac PCs differentiation Traditionally, the EB method has been widely used to examine the differentiation potential of PSCs in vitro. 17,30,31We set out to determine the credibility of the mESC reporter line during cardiac PCs differentiation using the EB method by monitoring fluorescent cells (Supplementary Fig. S2a).First, a drop of M10 medium containing 1000 mESCs was hung on the lid of a tissue culture dish.Five days later, the EBs formed and were plated into the gelatin-coated plate with WNT inhibitor IWR1-endo (days 5-7).On day 7, very few SHOX2 + ; HCN4 + cells emerged (Supplementary Fig. S2b).On day 9, about 1.6% of dual positive cells were determined by flow cytometry analysis, which lasted with a similar scale (1.31%-1.84%) in the whole differentiation process of the EB method (Fig. 1d and Supplementary Fig. S2c, d).Notably, on day 13, the expression pattern of Shox2:EGFP; Hcn4:mCherry EBs was similar to that of WT EBs.High expression of cardiac lineage marker Tnnt2, atrial markers Myh6 and Myl7, and ventricular marker Myl2 revealed that atrial-like cardiomyocytes (ALCMs) and ventricular-like cardiomyocytes (VLCMs) were the main cardiac cell types in EBs, whereas cardiac PC markers (Shox2 and Hcn4) 32,33 and cardiac PC-associated developmental gene Tbx18 were lowly expressed (Supplementary Fig. S2e-g).These data suggested by observing the fluorescent cells during differentiation, the cardiac PCs may be identified for further analysis.
The small quantity of cardiac PCs derived from the EB method encouraged us to develop a new, simple, and efficient method for cardiac PCs generation.A monolayer human PSC differentiation protocol with sequential WNT activation and inactivation has been applied to produce working cardiomyocytes with more than 90% efficiency. 34However, the efficiency of the monolayer cardiac PCs differentiation method from mouse PSCs is still very poor and needs to be optimized.We modified a previous EB method 30 by culturing Shox2:EGFP; Hcn4:mCherry mESCs in serum/LIF, which resembles a primed pluripotent state, as a monolayer (called M10 method) (Supplementary Fig. S3a).Cardiac PCs differentiation was induced by M10 medium supplemented with Vitamin C (Vc) for 3 days (days 0-3) and the cells were treated with WNT activator CHIR99021 (CHIR) for 24 h (days 2-3), then switched to serum-free medium for 2 days (days 3-5) and supplemented WNT inhibitor IWP2 for 2 days (days 5-7).On day 7, some HCN4 + cells emerged, but no SHOX2 + cells were observed (Fig. 1e and Supplementary Fig. S3b, c).After further differentiation, three subpopulations of cells with SHOX2 + ; HCN4 − , SHOX2 − ; HCN4 + , and SHOX2 + ; HCN4 + were identified in culture, all with extremely low ratios though (Fig. 1e and Supplementary Fig. S3b).The flow cytometry analysis detected only 0.52% and 0.41% SHOX2 + ; HCN4 + cells at days 11 and 13, respectively, during the differentiation period (Fig. 1e and Supplementary Fig. S3b, d), even less than the EB method (Supplementary Fig. S2c).Taken together, consistent with previous reports, 30 both the EB and M10 methods produced a low number of cardiac PCs, indicating the importance of developing a new cardiac PCs differentiation strategy.
The transition from naïve to formative pluripotency rendered efficient cardiac PCs differentiation of mESC reporter line As naïve PSCs are more homogenous than primed PSCs in terms of transcriptome and epigenome, 35 and FSCs are prone to lineage commitment, 22 we plated naïve 2i/LIF mESCs as seeds for cardiac PCs differentiation.The formative state (FS) transition was achieved by treating naïve mESCs with Activin A, tankyrase inhibitor XAV939, pan-retinoic acid receptor inverse agonist BMS493 and KnockOut serum replacement (KSR) for 3 days, evidenced by existence of pluripotency factor Oct4, increased expression of formative markers, Otx2, Oct6, Fgf5, and absence of naïve markers, Nanog, Tbx3, and Esrrb (Fig. 2a), which was similar to previous studies. 22,24The PSCs proliferated significantly during the 72-h FS transition.After this process, ubiquitous OTX2 expression in FSCs was observed (Supplementary Fig. S4a, b).Vc is included in media for both the EB and M10 methods (Supplementary Figs.S2a and 3a).Apart from its anti-oxidant effects, Vc has been shown as an important epigenetic modifier that plays a critical role in cell fate transition. 36We, therefore, added Vc after mesoderm commitment from the FS.Subsequently, after 2-day (days 3-5) Activin A, CHIR and Vc (called ACC) treatment, loss of pluripotency marker Oct4, upregulation of primitive streak markers Brachyury (T) and Foxa2, lateral cardiac mesodermal lineage-specific genes (Mesp1 and Flk1), and early cardiac lineage marker Nkx2-5 were detected, whereas the expression of endoderm markers, Sox17 and Hhex, was relatively low (Fig. 2b).We further examined the protein expression of T and MESP1 during mesoderm induction.The majority of cells on day 4 were T + and MESP1 + , which declined gradually with further induction, suggesting the efficient differentiation of cardiac mesodermal cells (Fig. 2c).Then, after continuous culture, 9.83 ± 5.81% SHOX2 + ; HCN4 + cells were generated on day 13 from ACC induction treated FS transition (FS + ACC), which was 32 folds more than the M10 method (Fig. 2d).
Small molecules targeting WNT, Activin/Nodal/TGF-β, and RA pathways further increased cardiac PCs differentiation efficiency of mESC reporter line As activation of RA and inhibition of WNT and Activin/Nodal/TGF-β signaling pathways have been reported to induce cardiac PCs differentiation, 14,[37][38][39] we first applied RA, A83-01 (TGF-β type I receptor inhibitor), and IWP2 individually during the differentiation of cardiac mesoderm into cardiac PCs from days 6 to 8. Five days after induction, the fluorescent cells were analyzed by flow cytometry, within the concentrations of tested small molecules, both RA and A83-01 increased SHOX2 − ; HCN4 + cells without effects on SHOX2 + ; HCN4 + cells.For IWP2, about an 11% increase of SHOX2 − ; HCN4 + cells were detected at high concentration (10 μM), however, there was about a 3% reduction of SHOX2 + ; HCN4 + cells as well (Supplementary Fig. S4c).
We then combined A83-01, IWP2, and RA together with Vc (called AIRC) for differentiation induction.Following FS conversion through FSK medium and combining with ACC induction and AIRC induction, we named the new cardiac PCs differentiation protocol as FSK method.
Although there was a slight increase in SHOX2 − ; HCN4 + cells (5%), an about 2.3-fold increase in SHOX2 + ; HCN4 + cells was achieved, reaching 17.57 ± 3.50%, compared with control (Fig. 2e), suggesting the synergetic effects of these signaling pathways in promoting cardiac PCs differentiation.A wide range of fluorescent cells was observed under the microscope on day 13 (Fig. 2f).We monitored the dynamic change of fluorescent cells by flow cytometry, which showed gradual increase of SHOX2 − ; HCN4 + and SHOX2 + ; HCN4 + cells from day 8 and reached apex on day 12, with 32.43 ± 3.03% and 22.1 ± 1.10% respectively, then SHOX2 − ; HCN4 + cells began to decrease, whereas SHOX2 + ; HCN4 + cells were relatively stable.On the other hand, less than 6% SHOX2 + ; HCN4 − cells were produced during the whole process (Fig. 2g and Supplementary Fig. S4d, e).
We thus established a new paradigm for cardiac PCs differentiation, the FSK method, which is composed of naïve to formative pluripotency transition, sequential mesodermal, and cardiac PCs induction (Fig. 2l).
We then performed cardiac PCs differentiation of WT E14TG2a mESCs by applying the FSK method.By immunofluorescence, compared to the M10 method, cells with high expression of cardiac PCs markers HCN4, SHOX2, and VSNL1 were widely distributed in day 13 differentiated cells accompanied by cTnT + and NKX2-5 + cells, further validating the robustness of this method (Supplementary Fig. S4h).These results implicated that our FSK method might induce a higher ratio of FSCs to differentiate into cardiac PCs.SHOX2 + ; HCN4 + cells displayed cardiac PC characteristics The fluorescence reporter cell line is beneficial for specific subtype purification and further functional characterization.To identify the cell types representing cardiac PCs, we sorted differentiated cells on day 13.The sorted cell morphologies displayed spider-like or spindle-like features (Supplementary Fig. S5a, b and Supplementary Video 1).Gene expression profile revealed that besides Shox2 and Hcn4, SHOX2 + ; HCN4 + cells enriched other cardiac PC markers, including Vsnl1, Hcn2, and Cacna2d2 (Fig. 3a and Supplementary Fig. S5c).The high expression of these ion channel and calcium-binding protein-encoding genes is the basis for the cardiac PCs' function. 33,40Tbx5 and Rgs6, as cardiac PC markers, were also enriched in SHOX2 + ; HCN4 + cells, suggesting the cardiac PC regulatory network was activated (Supplementary Fig. S5c).For SHOX2 + ; HCN4 − cells, the upregulation of Hcn2 may compensate for the relatively low Hcn4 expression (Fig. 3a).Although Cacna2d2, Tbx5, and Rgs6 were detected in SHOX2 − ; HCN4 + cells, lower expression of Vsnl1 32 and lack of Shox2 41 may compromise their cardiac PCs identity (Fig. 3a and Supplementary Fig. S5c).Meanwhile, in the context of the wide presence of early cardiac lineage markers (Myh6 and Nkx2-5) and cardiac cytoskeleton gene Tnnt2, we speculated that day 13 differentiated cells were similar to embryonic cardiac cells (Supplementary Fig. S5c).
Cardiac PCs have distinct electrophysiological properties.After cell sorting, we first performed a spontaneous action potential (sAP) assay and recorded diverse sAP forms representing three cell types, including cardiac PC, ALCM, and VLCM (Fig. 3b).In general, SHOX2 + ; HCN4 + cells displayed a fast spontaneous emission rate, slow maximum uplink speed, more positive maximum diastolic potential (MDP) and small action potential amplitude (APA) compared with other subpopulations (Fig. 3c and Supplementary Fig. S5d).The maximum uplink rate of SHOX2 + ; HCN4 + cells (38.53 ± 7.92 V/s) was 42.08%, 64.88%, and 84.44% of SHOX2 − ; HCN4 − (91.55 ± 14.23 V/s), SHOX2 + ; HCN4 − (59.39 ± 11.44 V/s) and SHOX2 − ; HCN4 + (45.63 ± 10.28 V/s) cells respectively (Fig. 3c).According to the results of sAP recording, on average, 45.31% of the SHOX2, HCN4 single and double positive cells displayed the functional sAP of cardiac PCs (Fig. 3d), which meant a large portion of cardiac PCs can be obtained from the FSK method.55.00% of SHOX2 + ; HCN4 + cells showed cardiac PC sAP, which was the highest among the four sorted populations.A similar ratio of ALCM sAP exists in two subpopulations (48.15% in SHOX2 + ; HCN4 − cells and 47.06% in SHOX2 − ; HCN4 + cells) while 36.36% of SHOX2 − ; HCN4 − cells were considered as VLCMs (Fig. 3d).To further characterize the electrophysiological properties of SHOX2 + ; HCN4 + cells, the I f current was also recorded in SHOX2 + ; HCN4 + cells which were sensitive to ivabradine (Fig. 3e).We further detected the effect of cAMP on I f current in SHOX2 + ; HCN4 + cells.Similar to SAN cells, I f currents showed an increase in response to cAMP activity (Supplementary Fig. S5e, f).Meanwhile, SHOX2 + ; HCN4 + cells could maintain cardiac PC characteristics at least up to day 30 with similar maximum uplink rate and ivabradine-sensitive I f currents (Supplementary Fig. S5g-i).In all, these data indicated that SHOX2 + ; HCN4 + cells were more similar to cardiac PCs.
Next, we examined the pacemaker potential of SHOX2 + ; HCN4 + cell aggregates in vitro.The cells were plated into a lowattachment U-bottom 96-well plate and centrifuged to make aggregates (Supplementary Fig. S5j).The aggregates were then seeded onto neonatal mouse ventricular cardiomyocytes (NMVMs), and the pacing activity of cardiac PC aggregates were assessed by multiple microelectrode array (MEA) assays.Generally, monolayer NMVMs alone have multiple initiation sites.In the presence of SHOX2 + ; HCN4 + cardiac PCs aggregates, a dominant pacemaker stably paced NMVMs, increasing the pace from 155 beats per minute (b.p.m.) to 175 b.p.m., and generated more rhythmic beats; higher FPD max (FPD, field potential duration) and larger spike amplitude were detected as well (Fig. 3f, g).These results indicated that SHOX2 + ; HCN4 + cells functioned as cardiac BP in vitro.
Finally, we tested the pacemaker activity of differentiated cells in vivo (Supplementary Fig. S5k).We injected 1-2 × 10 6 differentiated cells into the apex of the rat heart with daily immunosuppressive treatment.Optical mapping revealed that the implanted cells paced the rat's left ventricle for a short period during drug-induced short-term atrioventricular conduction block (Fig. 3h).Subsequent paraffin section confirmed the presence of implanted cells in the heart by co-localization of SHOX2 and HCN4:mCherry, which may form a connection with host ventricular cardiomyocytes via CX43 (Supplementary Fig. S5l, m).
Cardiac PC cluster was identified by scRNA-seq during in vitro cardiac PCs differentiation To unravel the dynamic process of cardiac PCs differentiation induced by the FSK method, we utilized a droplet-mediated scRNA-seq platform (10× Genomics Chromium) to capture cells at different time points (days 5, 8, and 13), corresponding to the completion of mesoderm, cardiac PCs induction, and later stage maintenance respectively (Fig. 4a).We sequenced 81,599 cells at an average of about 50,000 reads per cell (Supplementary Table 1).After quality control, cells with aberrant gene detections (< 500 or > 6000 genes) were removed.Especially, we filtered cells with high mitochondrial gene coverage (> 10%) in day 5 samples Fig. 3 SHOX2 + ; HCN4 + cells displayed cardiac PC characteristics.a qRT-PCR analysis of the expression of cardiac lineage and cardiac PC marker genes in sorted cells.Relative to Gapdh expression, normalized to SHOX2 − ; HCN4 − cells.mCherry and EGFP were closely correlated to Hcn4 and Shox2 expression, respectively.Cardiac PC-related TF marker: Shox2; cardiac PC-related ion channels: Hcn4, Hcn2, and Cacna2d2; cardiac PC enriched gene: Vsnl1.Red, SHOX2 + ; HCN4 + cells; yellow, SHOX2 + ; HCN4 − cells; deep blue, SHOX2 − ; HCN4 + cells; light blue, SHOX2 − ; HCN4 − cells.Data were presented as means ± SEM from technical triplicates (n = 3).b Representative sAP recording of single cells differentiated from Shox2:EGFP; Hcn4:mCherry mESCs via the FSK method.c The maximum uplink speed of four different subpopulations.Red, SHOX2 + ; HCN4 + cells (n = 20); yellow, SHOX2 + ; HCN4 − cells (n = 27); deep blue, SHOX2 − ; HCN4 + cells (n = 17); light blue, SHOX2 − ; HCN4 − cells (n = 22).Data were presented as means ± SEM, n indicated the number of recorded cells.d The proportion of different cell types according to AP forms in the sorted subpopulations.Red, cardiac PC; yellow, ALCM; light blue, VLCM.e The I f current recording of SHOX2 + ; HCN4 + cells.Baselines were recorded before adding 1 µM ivabradine into Tyrode's solution with 1 mM Ba 2+ .f Activation maps of electrical signal propagation revealed that impulses generated from cardiac PC aggregates.Colormap indicated the origin of impulses generated from the electrode 11.Cardiac PC aggregate was laid next to electrode 11 and propagated electrical signal to pace NMVMs.Red arrows represented the directions of impulse propagation.g Representative raw traces showed the addition of cardiac PC aggregate increased the beating rate.Representative raw traces of the monolayer NMVMs with cardiac PC aggregate were recorded at electrodes 11, 33, and 44 with the speed of 175 b.p.m. h Differentiated cells from the FSK method acted as an ectopic pacemaker in rat atrioventricular conduction block model.Yellow stars indicated the pacemaking sites.Cells were sorted and collected on day 13 differentiation (a-e).sAP recording and I f current recording were performed on HEKA EPC-10 amplifier (b-e).Comparisons between multiple groups were performed with a one-way ANOVA test (a) and Kruskal-Wallis rank sum test (c).Statistical significance was indicated as follows: ns not significant; P < 0.05 (*); P < 0.01 (**); P < 0.001 (***); P < 0.0001 (****) because quality control of mitochondrial RNA (mtRNA) might introduce a bias that particularly discriminates cardiac PCs (days 8 and 13). 42In this manner, we sequenced 66,681 cells at three time points, with 22,527, 11,236, and 28,918 in day 5, 8, and 13 samples, respectively (Fig. 4b).Unsupervised clustering of the scRNA-seq data identified 23 distinct clusters (cluster 0-22, Supplementary Fig. S6a, b), which were annotated into 16 major cell types based on cell-specific markers (Fig. 4c, d As expected, the majority of cells on day 5 were mesodermal cells (Fig. 4e and Supplementary Fig. S6c).The clusters of Pacemaker progenitor cells and Cardiac PCs already emerged on day 8, consistent with our experimental observations (Fig. 2g).The number of mesoderm and pacemaker progenitor cells was significantly reduced on day 13 (Fig. 4e and Supplementary Fig. S6c), reflecting the sufficient differentiation from mesoderm to cardiac PC lineage.We next analyzed the expression of two cardiac PCs markers (Shox2 and Hcn4).Although both Shox2 and Hcn4 were detected on day 8, the expression of both markers was higher on day 13 (Fig. 4f), consistent with the flow cytometry data, indicating that more cardiac PCs were generated and some cardiac PCs might be more mature on day 13.Besides, genes upregulated in the Cardiac PC cluster were enriched in the function of cardiomyocyte differentiation regulation, muscle contraction, ion homeostasis, and circadian entrainment, etc. (Fig. 4g).The Cardiac PC cluster not only highly expressed Shox2 and Hcn4 but also expressed abundant SAN and cardiac PC markers mentioned in previous research, such as Vsnl1, Bmp4, Cacna2d2, Igfbp5, Lbh, and Dact1 (Fig. 4h and Supplementary Fig. S6d). 32,43Notably, the Shox2 + ; Hcn4 + cells were mainly in the Cardiac PC cluster, while the Hcn4 + cells were not restricted to this cluster (Fig. 4h), suggesting that Hcn4 single positive might import a false-positive definition of cardiac PCs.To compare our differentiated cardiac PCs with the in vivo SAN cells, we integrated our data with the scRNA-seq data of E13.5 mouse SAN cells (GSE130461, Supplementary Fig. S6e). 44As indicated by the expression of Shox2 and Hcn4 markers, SAN cells in mice were highly overlapped with our differentiated cardiac PCs in Uniform Manifold Approximation and Projection (UMAP), illustrating that our in vitro differentiated cardiac PCs obtained from the FSK method were similar to in vivo PCs.Together, these results showed that dual positive markers (Shox2 and Hcn4) could define cardiac PCs more precisely, and our in vitro cardiac PCs differentiation protocol, the FSK method, guided a credible pacemaker lineage.

Cardiac PCs differentiation trajectory was successfully charted
To explore the refinements in cardiac PCs differentiation over time, we employed pseudo-time analysis to construct the differentiation trajectory with our scRNA-seq data (Fig. 5a).The trajectory depicted three different states, bifurcating from the main pre-branch (State 1) into two branches representing a successful branch to cardiac PCs (State 3) and a branch to cardiac fibroblasts (State 2) (Fig. 5a, b).Gsc, a mesoderm marker, was highly expressed in the initial stage of State 1 (Fig. 5c).The expression of Isl1 and Tbx18 at the crossroad of cardiac PCs differentiation trajectory suggested that the progenitor cells with posterior heart field (PHF) characters might contribute to cardiac PCs generation. 45Two cardiac PC markers (Shox2 and Hcn4) indicated the successful branch to cardiac PCs in State 3, while the high expression of Col1a1 represented another branch to cardiac fibroblasts (Fig. 5c).
To comprehensively view the differentiation trajectory, we analyzed 2829 differentially expressed genes (DEGs, with qvalue < 0.01) and observed six gene expression clusters with different patterns (Fig. 5d).Genes in expression cluster 1 (C1) (a gene cluster referred to the mesoderm stage), enriched in ribosome biogenesis, were gradually downregulated from prebranch (State 1), showing the translation activity was declining during the differentiation process (Fig. 5d).Conversely, the proliferation and sterol biosynthesis processes were upregulated (C2 and C3), especially in the middle stage of the differentiation trajectory (Fig. 5d), which might be due to the existence of progenitor cells, such as CPCs and pacemaker progenitor cells.Genes that increased along the State 3 trajectory were involved in the functions of NADH dehydrogenase complex assembly and muscle contraction (Fig. 5d), consistent with the nature of cardiac PCs. 46,47Notably, the failed branch (State 2) was highly enriched in extracellular matrix (ECM) organization (Fig. 5d, e), suggesting that the trajectory to cardiac fibroblasts might be a branch during our cardiac PCs differentiation.Since TFs are critical in the specification and differentiation of cardiac PCs, 48 we further analyzed the expression of TFs in six gene expression clusters (Fig. 5d).TFs expressed in the mesoderm, such as Gsc, Pitx2, and Hes1, were highly expressed in C1.Several progenitor-related TFs (Sox11, Ets1, Gata5, and Gata6) were enriched in C2 and C3.Tcf4, Zeb2, Zfp503, and Id2 were also highly expressed in C3, indicating their potential roles in cardiac PCs differentiation.C5 represented the state of cardiac PCs and consisted of abundant TFs expressed in PCs, such as Mef2c, Shox2, and Tbx5.Moreover, Id1, Jun, Fos, Klf6, Gata4, and Zeb1 in C6 were closely related to the process of ECM, reflecting its trajectory to the cardiac fibroblasts.Figure 5e shows the representative TFs during the cardiac PCs differentiation.In summary, these data depicted the trajectory of cardiac PCs differentiation and revealed potential critical TFs for PCs development.
ZFP503 functioned effectively in the cardiac PCs differentiation While reviewing the hit gene list from the single-cell transcriptome, a zinc finger TF, Zfp503, came to our attention.Zfp503 barely expressed on day 5 while significantly increased on days 8 and 13 (Fig. 6a, b).The dynamic change of Zfp503 expression in cardiac PCs differentiation was consistent with the pseudo-time analysis (Figs. 5e and 6c), which displayed that Zfp503 expressed in the early period of the cardiac PCs developmental process, from SHF, CPCs, pacemaker progenitor cells to cardiac PCs (Fig. 6d).As an RA-related gene, Zfp503 significantly enriched in the Cardiac PC cluster (Fig. 6d).Also, from scRNA-seq of E13.5 mouse SANs, Zfp503 was enriched in primary PCs, implying its potential role in SAN development both in vivo and in vitro.High level of Zfp503 in sorted SHOX2 + ; HCN4 + cells was validated by qRT-PCR and immunofluorescence staining (Fig. 6e and Supplementary Fig. S7a).Concomitantly, ZFP503 co-expressed with HCN4 in Fig. 4 Cardiac PC cluster was identified by scRNA-seq during in vitro cardiac PCs differentiation.a Schematic representation of the timeline of sample collection for scRNA-seq (days 5, 8, and 13) (n = 2).b UMAP dimensionality reduction plot of samples on days 5, 8, and 13. 22,527, 11,236, and 28,918 cells were sequenced from days 5, 8, and 13 samples, respectively.c UMAP and clustering of scRNA-seq data identified 16 distinct clusters (VLCM, TC, SHF, Cardiac PC, Primitive streak-derived tissues, Paraxial mesoderm, Pacemaker progenitor cell, Nascent mesoderm, Mesendoderm/epiblast, Immature cardiomyocyte, Fibroblast, Endothelial, Endoderm, CPC, Cardiac fibroblast, and ALCM).Each dot represented an individual cell colored by an annotated cluster.TC transitional cell; SHF second heart field.d Dot plot of cell-specific markers defining each cluster.e UMAP dimensionality reduction plots divided by samples on days 5, 8, and 13.Different colors represent scRNA-seq clusters defined in (c, d).f Violin plots of the normalized unique molecular identifier (UMI) counts for two cardiac PC markers (Shox2 and Hcn4) in days 5, 8, and 13 samples.Each dot in the plot represents one cell.The expression level means the log normalized data.Kruskal-Wallis rank sum test: P < 0.001 (***).g GO (gene ontology) term analysis for the genes upregulated in the Cardiac PC cluster.GO terms with P < 0.05 was considered significant.h UMAPs depicting the expression density of typical PC markers genes (Hcn4, Shox2, Bmp4, and Vsnl1) upregulated in the Cardiac PC cluster  6f and Supplementary Fig. S7b).
Previous findings defined Zfp503 as a marker of neural progenitor cells and an RA-activated gene. 25,49The putative RA response element (RARE) is located 1600 bp downstream of the Zfp503 transcriptional start site (TSS). 25Luciferase reporter assays further confirmed that RA activated the transcriptional activity of Zfp503 (Supplementary Fig. S7c).While the RA signaling pathway is important in cardiac PCs' fate determination, the relationship and mechanism of Zfp503 and cardiac PCs differentiation need to be further elucidated.Not surprisingly, a low dose of RA (0.25 μM) activated Zfp503 expression, whereas AIRC treatment (containing 0.25 μM RA) promoted a similar degree or slight increase of Zfp503 expression (Fig. 6g).Zfp503 was activated by RA in a dosagedependent manner, which was repressed by BMS493 (Fig. 6g and Supplementary Fig. S7d).These data confirmed Zfp503 as an RA target gene in the cardiac PCs differentiation.
We speculated that Zfp503 respond to the RA signaling pathway and functions in the process of cardiac PCs development.Zfp503 increased with RA treatment and functioned in the cardiac PCs developmental process and promoted the generation of SHOX2 + ; HCN4 + cells.RA boosted the effect of A83-01, IWP2, and Vc (called AIC) treatment and upregulated the efficiency of SHOX2 + ; HCN4 + cell differentiation (Supplementary Fig. S7e).

DISCUSSION
High-quality seeds are paramount for qualitative production.So far, most mouse PSCs differentiation methods start with cells cultured in serum/LIF, which are heterogeneous with cells biased to differentiation lineages. 50,51While various approaches have been attempted to optimize the signaling pathways to direct the differentiation routes, the quality of 'seeds', i.e., PSCs, was largely overlooked due to the limited culture media for PSCs.Recently, it has been proposed and verified along the embryo development, upon the exit from naïve pluripotency, the epiblast transiently goes through a FS to acquire differentiation potential before committing to a certain lineage.It has been demonstrated that FSCs could be induced to differentiate efficiently. 22,52Based on this progress in stem cell research, we established a new mouse cardiac PCs differentiation method by transitioning naïve PSCs to formative pluripotency for subsequent directed cardiac PCs differentiation, which produced SHOX2 + ; HCN4 + cells with typical PC characteristics, showcasing the applicability of our strategy.
4][55][56] Manipulating WNT signaling has been widely used in cardiac lineage differentiation from PSCs.Due to the difference in the application of WNT activation and suppression during various cardiac PCs differentiation methods, the critical role of WNT signaling in cardiac PCs differentiation is still not clearly defined. 31,34,38,57Yechikov et al. showed that inhibition of Activin/Nodal/TGF-β signaling pathway together with WNT pathway improved the efficiency of cardiac PCs differentiation. 58RA signaling pathway was also demonstrated to be critical in promoting atrial and SAN-like fate while inhibiting ventricularlike differentiation. 31,55,59Although Protze et al. showed that RA increased the level of cardiac PC markers but did not affect cardiac PCs differentiation efficiency, the use of RA in conjunction with BMP during mesoderm induction promoted cardiac PCs differentiation.Many other studies have also demonstrated that RA promoted the differentiation of PCs when combined with other signaling modulators. 55,59Therefore, we decided to include RA together with A83-01 and IWP2 in our protocol.
Interestingly, in our observation, the effects of RA, A83-01, and IWP2 alone were not significant in obtaining SHOX2 + ; HCN4 + cells.However, when combined, the number of SHOX2 + ; HCN4 + cells increased significantly, indicating the complex and delicate regulation of signaling pathways in cardiac PCs differentiation.On the other hand, the crosstalk among these pathways may direct the differentiation route of cardiac PCs 59 , whereas manipulating a single pathway was not significant.We, therefore, optimized the time and strength of either activating or repressing these pathways during cardiac PCs differentiation with small molecule combinations, which further increased the dual positive cells to more than 20%.By applying more specific and potent chemicals targeting the relevant pathway/component at the right differentiation window, more efficient cardiac PCs differentiation may be realized.
The conservation of gene networks in SAN of different species including mouse and human has been characterized. 32,33However, physiologically, the heart rates for mice and humans are about 670 and 70 b.p.m., which is significantly different.The spontaneous rate is 4.6-times faster, and the action potential (AP) is 2.4-times shorter, while the I f current is about 5-fold higher in mouse than in human, which may be caused by higher expression of ion channels such as HCN1, 2, 4, Cav3.1, RyR2 and SERCA2 in mouse SAN. 62Several laboratories have done pioneering work on human cardiac PCs (hPCs) differentiation that also laid a solid foundation for our current research.By modulating BMP, Activin/ Nodal/TGF-β, and RA signaling pathways, Protze et al. induced human ESC (hESC) (HES3-NKX2-5 gfp/w ) into SAN-like pacemaker cells (SANLPCs), characterized by positive for pan-cardiomyocyte surface marker SIRPA and mesenchyme marker CD90, while negative for NKX2-5:GFP. 31More recently, Han et al. established SHOX2:GFP; MYH6:mCherry hESC line to produce hPCs. 37,61Using nine substances targeting singling pathways and epigenetic modifiers, about 46.6% SHOX2:GFP + cells were produced, though only SHOX2 is the cardiac PC-specific expression gene in this reporter system. 37Wakimizu et al. generated HCN4-EGFP transgenic, SHOX2-mCherry knock-in reporter human iPSC line, with different hPCs differentiation protocol, the dual positive cells ranged from 34%-44%.The gene editing strategies applied in this study may result in an N-terminal fusion protein or disruption of the coding sequence of HCN4 and SHOX2. 63In the future, more cardiac PC-specific expression genes or gene combinations are needed for PCs purification and characterization.
sAP of cardiac PCs is the key to the development of cardiac BP technology.In this study, AP recording revealed that 55.00% SHOX2 + ; HCN4 + cells had typical cardiac PC-like APs and characteristic I f current, which could increase by cAMP (Supplementary Fig. S5e, f) and inhibit by ivabradine (Fig. 3e and Supplementary Fig. S5i).Interestingly, in SHOX2 + ; HCN4 − and SHOX2 − ; HCN4 + subpopulations, 44.44% and 35.29% cells displayed cardiac PC-like APs, respectively (Fig. 3d).It is worth noting that there was still a 22.73% population of cardiac PCs in the SHOX2 − ; HCN4 − population.This might be attributed to our stringent cell sorting gates to obtain a relatively homogenous cell population.After cell sorting, low expression of mCherry, Hcn4, and Vsnl1 was detected in SHOX2 − ; HCN4 − cells, which may express low levels of protein to enable some cells to display cardiac PC-like APs.These data suggested that using one cardiac PC marker alone is not reliable in PCs identification; on the other hand, by tracing the cell fate transition among these single and dual positive populations, we may further improve the cardiac PCs differentiation efficiency.In addition, we demonstrated that by MEA, cardiac PC aggregates acted as a dominant pacemaker to pace primary NMVMs.More critically, after cardiac implantation of cardiac PCs, they were capable of ectopically pacing the ventricle of rat heart ex vivo.
scRNA-seq facilitates us in understanding developmental trajectories and key factors determining cardiac PCs differentiation.The day 5 dataset revealed that most cells were mesendoderm progenitors, verifying the efficient commitment of FSCs after activation of WNT and Activin/Nodal/TGF-β signaling pathway. 22n subsequent small molecule cocktail-directed cardiac PCs differentiation, the trajectory showed that a portion of cardiac PCs was developed from Isl1 and Tbx18 expressing cells, indicating the in vitro development process of cardiac PCs may originate from progenitor cells with PHF properties, mirroring the development of bona fide cardiac PCs in vivo (Fig. 5c). 64,65portantly, we validated Zfp503, a downstream target of RA signaling pathway, 49 functioning in the development of cardiac PCs.It has been documented that Zfp503 is involved in regulating early neuromesodermal progenitors and later brain and limb development. 25,66,67Zfp503 is a potential player mediating the crosstalk of BMP, WNT, and RA signaling pathways. 67In our scRNAseq dataset, Zfp503 was expressed from SHF progenitor cells to cardiac PCs, which highly matched the developmental process of cardiac PCs.Moreover, Zfp503 null mutation resulted in low efficiency of SHOX2 + ; HCN4 + cardiac PCs differentiation.RA appropriately induced Zfp503 expression and promoted the differentiation of SHOX2 + ; HCN4 + cardiac PCs, consistent with previous reports that the RA signaling pathway promoted the development of cardiac PCs. 39,59he top DEGs of each cluster also indicated the co-existence of different cardiac progenitors.These results, on the one hand, suggested the possibility of efficiently obtaining atrial and ventricular cardiomyocytes with a modified FSK method.On the other hand, it will further direct us to improve the cardiac PCs differentiation efficiency.For example, Nkx2-5 was highly expressed in the first heart field (FHF) and SHF clusters.It has been demonstrated the antagonism between Nkx2-5 and Shox2 in cardiac progenitors determined their commitment either to cardiomyocytes or to cardiac PCs. 9 Therefore, by manipulating the cell fate transition, the yield of cardiac PCs may increase further.
In summary, we have established a simple and efficient platform for the generation of mouse cardiac PCs from PSCs.With scRNA-seq, we mapped the in vitro development trajectory of cardiac PCs and verified Zfp503 as a key factor in guiding cardiac PCs differentiation.These findings will not only provide a high-quality cell source for cardiac BP but enhance our understanding of cardiac PCs development as well.By establishing specific reporter PSC lines, following the developing trajectory, and modulating the relevant signal pathways, our integrated cardiac PCs differentiation strategy may be extended to other cell types that are urgently needed in regenerative medicine.

Ethics statement
All the animal experiments were approved by the Animal Care and Use Committee at the School of Medicine, Tongji University (No. TJBB00921701).Neonatal mice were euthanized by decapitation before harvesting hearts for cardiomyocyte isolation.For in vivo pilot transplantation, rats were intubated for mechanical ventilation followed by anesthesia with 4% isoflurane.For the isolation of adult hearts, rats were euthanized by cervical dislocation after anesthesia with 4% isoflurane.All procedures were performed in accordance with the Guide for the Care and Use of Laboratory Animals made by the U.S. National Institutes of Health.Animals Neonatal C57/BL6J mice (1-3 day-old) and Sprague Dawley rats (male, 150-200 g, 6-week-old) were used in this study.All animals were purchased from Shanghai Sippe-Bk Lab Animal Co., Ltd., China.
mCherry was first introduced into the Hcn4 locus.After another 4 to 6 days, individual colonies were picked, expanded, and characterized by PCR amplification and subsequent sequencing for correct targeting.The same process was performed to introduce 2A-EGFP into the Shox2 locus in Hcn4:mCherry mESCs.
Embryoid body (EB) method mESCs cultured in serum/LIF were plated as hanging drops on the lid of a tissue culture dish with M10 medium (KnockOut DMEM plus 10% fetal bovine serum) supplemented with 0.5 mM Vc.Every drop contains 1000 cells in a volume of 20-30 µL.At day 5, EBs were collected and plated into a gelatin-coated 6-well plate.5 µM IWR1-endo and 0.25 mM Vc were added into the M10 medium for 2 days to induce cardiac differentiation.The medium was switched to maintain medium (MEM supplemented with Insulin-Transferrin-Selenium-Sodium Pyruvate, penicillin streptomycin glutamine, and β-mercaptoethanol) for short-term cell culture (7-14 days) before analysis.

M10 method
In total, 2 × 10 4 mESCs cultured in serum/LIF were plated into one well of gelatin-coated 6-well plate overnight.The next day, the medium was switched to M10 containing 0.5 mM Vc. Two days later, CHIR (3 μM) was added to activate the WNT signaling pathway for 24 h.RPMI1640 supplemented with B27 minus insulin and 0.5 mM Vc (RBC-i) was used to differentiate the cells for 2 days.Then, cells were treated with 5 μM WNT inhibitor IWP2 for 48 h to promote cardiac lineage differentiation.The cells were maintained in RPMI1640 with B27 and Vc (RBC) till further analysis.

FSK method
In total, 5 × 10 5 naïve mESCs were plated into one well of fibronectin (10 µg/mL)-coated 6-well plate.The next day, the medium was switched to FSK medium containing N2B27 supplemented with 2 µM BMS493, 2 µM XAV939, 3-6 ng/ mL Activin A, and 1% KSR.After 72 h, the cells were gently dissociated and replated at 5 × 10 5 cells/well into fibronectincoated 6-well plate, induced with mesoderm derivation medium, including N2B27 containing 20 ng/mL Activin A, 0.5 mM Vc and 3 µM CHIR (ACC) for 48 h.Subsequently, the medium was changed to RBC-i.After 24 h, cells were treated with 5 µM A83-01, 5 µM IWP2, 0.5 mM Vc, and 0.25 µM RA (AIRC) for 2 days with everyday medium change.After that, the medium was changed back to RBC-i for 48 h with everyday medium change.In order to promote cell proliferation, cells were cultured in RBC from day 10 onward, and then harvested on days 13 or 30 for subsequent analysis.

Statistical analysis
All data represented three biological or technical replicates or more unless otherwise indicated in the figure legends.Statistical analyses were performed by GraphPad Prism Version 8.2.1(279) or R software (version 4.1.1).All data were presented as mean ± SEM.Error bars were indicated in the figure legends.All statistical analyses were conducted using a two-tailed student t-test, Mann-Whitney test, one-way ANOVA test, or Kruskal-Wallis rank sum test, where appropriate.P value < 0.05 was considered statistically significant.Statistical significance was indicated as follows: ns not significant; P < 0.05 (*); P < 0.01 (**); P < 0.001 (***); P < 0.0001 (****).

Fig. 5
Fig. 5 Cardiac PCs differentiation trajectory was successfully charted.a Unsupervised transcriptional trajectory of cardiac PCs differentiation, colored by pseudo-time and cell states (States 1-3).b Trajectory reconstruction of cardiac PCs differentiation revealed three branches: pre-branch (State 1), failed branch (State 2), and successful branch (State 3).Colors represented different cell clusters.c The expression of marker genes of three branches.Mesoderm markers: Gsc; cardiac PC markers: Hcn4, Shox2, Tbx18, and Isl1; cardiac fibroblast marker: Col1a1.d The differentially expressed genes along the pseudo-time were divided into six clusters (C1-C6) showing different patterns.The top GO terms of each cluster were shown.Highly expressed TFs of each gene expression cluster were shown.e The expression of representative genes (Nkx2-5, Otx2, Klf6, Nfat5, Zeb1, Zfp503, Shox2, and Tbx20) in the in the pseudo-time trajectory.The expression level means the log normalized data